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ABSTRACT 

The  oceanic  turbulent  boundary  layer  is  a  critical  region  to 
understand  for  oceanic  and  atmospheric  prediction.  This  thesis 
answers  two  fundamental  questions: (i)  what  is  the  response  of  the 
ocean  mixed  layer  system  to  transient  forcing  at  the  air  sea 
surface?  (ii)  what  is  the  necessary  time  and  space  resolution  in  an 
ocean  mixed  layer  model  to  resolve  important  transient  responses? 

Beginning  with  replication  of  de  Szoeke  and  Rhines'  (1976) 
work,  additional  physical  processes  were  added  to  include  more 
realistic  viscous  dissipation  and  anisotropy  in  the  three- 
dimensional  turbulent  kinetic  energy  (TKE)  budget.  These  refinements 
resulted  in  modification  of  de  Szoeke  and  Rhines'  findings. 
Firstly,  TKE  unsteadiness  is  important  for  a  minimum  of  105  seconds. 
Secondly,  viscous  dissipation  should  not  be  approximated  as  simply 
proportional  to  shear  production.  Thirdly,  entrainment  shear 
production  remains  significant  for  a  minimum  of  one  pendulum-day. 

The  required  temporal  model  resolution  is  dependent  on  the 
phenomena  to  be  studied.  This  study  focused  on  the  diurnal, 
synoptic,  and  annual  cycles,  which  the  one-hour  time  step  of  the 
Naval  Postgraduate  School  model  adequately  resolves.  The  study  of 
spatial  resolution  showed  unexpectedly  that  model  skill  was 
comparable  for  1  m,  10  m  and  even  20  m  vertical  grid  spacing. 
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I .   INTRODUCTION 

A.  STATEMENT  OF  PROBLEM 

The  oceanic  turbulent  boundary  layer  or  mixed  layer  is  a 
critical  region  to  understand  when  studying  the  atmosphere  or 
ocean.  The  ocean  is  the  primary  source  of  moisture  fluxed 
into  the  atmosphere,  it  is  a  significant  source  of  heat  fluxed 
into  the  atmosphere  in  many  regions,  and  it  is  a  moderator  of 
atmospheric  extremes.  To  accurately  model  the  atmosphere  or 
ocean,  a  realistic  representation  of  the  turbulent  mixed  layer 
is  reguired.  This  thesis  will  address  two  fundamental 
guestions: (i)  what  is  the  response  of  the  ocean  mixed  layer 
system  to  transient  forcing  at  the  air  sea  surface?;  (ii)  what 
is  the  necessary  time  and  space  resolution  in  an  ocean  mixed 
layer  model  to  resolve  important  transient  responses?  With 
the  answers  to  these  guestions,  air-ocean  models  will  be  a 
step  closer  to  being  realistically  coupled.  The  result  will 
be  better  predictions  for  both  atmosphere  and  ocean. 

B .  METHOD 

To  accomplish  this  task,  the  work  of  de  Szoeke  and 
Rhines(1976)  on  wind  driven  mixing  was  replicated  and  then 
generalized  to  include  surface  heating  and  cooling.  Their 
study  in  "Asymptotic  regimes  in  mixed  layer  deepening"  did  not 
include  physical  processes  associated  with  surface  heat  flux 
and  anisotropy  in  the  three-dimensional  turbulent  kinetic 


energy (TKE)  budget.  With  the  addition  of  these  critical 
physical  processes,  a  more  generalized  study  of  the  short  time 
scale  dependence  of  the  mixed  layer  is  possible.  The  study 
can  be  extended  past  the  single  day  scale  of  de  Szoeke  and 
Rhines  to  the  synoptic  and  annual  time  scales  by  shifting  the 
method  of  solution  from  a  Runge-Kutta  numerical  scheme  to  a 
FORTRAN-based  gridded  model.  With  the  question  of  the  minimum 
time  scale  resolved,  new  assumptions  about  mixed  layer  physics 
can  be  made  which  allow  testing  of  spatial  resolution.  The 
Naval  Postgraduate  School (NPS)  mixed  layer  model  developed  by 
Garwood  (1977)  was  modified  to  test  the  effect  on  model  skill 
of  varying  the  vertical  grid  size  using  observed 
meteorological  forcing  and  changing  thermal  structure  at 
weather  station  Papa  (50°  N,  145°  W)  . 


II.  THEORY,  MODEL  EQUATIONS,  AND  HYPOTHETICAL  SOLUTIONS 

A.  AIR  OCEAN  COUPLING 

The  surface  layer  of  the  ocean  responds  to  variations  in 
wind  stress,  heat  flux,  moisture  flux,  and  several  other  less 
significant  factors.  These  factors  have  been  examined  in 
depth  since  Ekman' s  (1905) ground  laying  work. 

B.  REPLICATION  OF  DE  SZOEKE  AND  RHINES  (1976) 

Initial  modeling  efforts  for  this  study  focused  on 
reproducing  and  then  building  on  the  efforts  of  de  Szoeke  and 
Rhines  (1976)  as  discussed  in  their  paper  "Asymptotic  regimes 
in  mixed-layer  deepening."  Niiler  (1975)  proposed  the  basic 
terms  in  the  turbulent  kinetic  energy  (TKE)  equation  studied 
by  de  Szoeke  and  Rhines: 

(i^2i24_u4_2_  (1_cosft)  +rufh2)  —  =2mnu?h  (1) 

2  f2  °       dt  ° 

A  B  CD 

Term  A  represents  buoyant  damping  of  TKE  by  entrainment,  term 
B  the  shear  production  of  TKE  by  entrainment,  term  C  the  spin 
up  or  storage  of  TKE,  and  term  D  the  near-surface  wind  shear 
production  of  TKE  minus  viscous  dissipation.  Solar  heating 
and  surface  heat  flux  were  not  included.  Term  C,  the  energy 
required  to  spin  up  the  turbulence  intensity,  was  not 
originally  included  by  Niiler. 


To  mathematically  define  the  components  of  the  de  Szoeke 
and  Rhines  model,  the  friction  velocity  or  u.  must  first  be 
defined.  The  friction  velocity  (equation  2)  is  defined  as 
the  square  root  of  the  surface  wind  stress  divided  by  the 
density  of  water. 


u  = 

'  N 


(0) 


The  Brunt-Vaisala  or  buoyancy  frequency  (N)  is  a  measure  of 
vertical  stability.  A  value  of  2II/600  sec-1  was  used  for 
initial  model  testing  and  in  the  hypothetical  cases.  The 
mixed  layer  depth  is  h.  The  Coriolis  parameter  (f)  is  defined 
as  2Q,  where  Q2  is  the  vertical  component  of  planetary 
rotation.  A  nominal  mid-latitude  value  of  lxlO-4  sec-1  was 
used  for  f .  In  equation  (1)  C0  and  m0  are  dimensionless  tuning 
coefficients  for  TKE  spin  up  and  net  wind-shear  production, 
respectively. 

Buoyant  damping  of  TKE  by  entrainment,  term  A,  is  the 
rate  of  decrease  in  turbulence  due  to  the  lifting  upward  of 
denser  water  originating  below  the  mixed  layer.  The  mixed 
layer  is  made  more  stable  by  the  entrainment  of  colder  denser 
water,  with  a  commensurate  increase  in  potential  energy. 
Shear  production  due  to  entrainment,  term  B,  is  the  conversion 
of  mean  kinetic  energy  to  turbulent  kinetic  energy  due  to 


shear  across  the  base  of  the  mixed  layer.  This  term  also 
results  in  a  net  gain  in  turbulence.  Turbulence  spin  up,  term 
C,  is  the  TKE  build  up  required  before  entrainment  can  occur. 
The  wind  stirring  minus  viscous  dissipation,  term  D, 
represents  the  net  input  of  energy  into  the  mixed  layer  due 
to  wind  stress,  i0,  in  equation  (2)  .  In  de  Szoeke  and  Rhines' 
model,  viscous  dissipation  is  also  proportional  to  u*3  because 
surface  heating  and  cooling  are  neglected  as  sources  or  sinks 
of  TKE,  and  dissipation  of  TKE  produced  by  entrainment  shear 
production  is  not  accounted  for. 

As  a  first  step  in  demonstrating  and  understanding  the 
physics  that  drive  mixed  layer  deepening,  de  Szoeke  and 
Rhines'  results  were  replicated.  Equation  (1)  was  solved  for 
the  rate  of  mixed  layer  deepening,  We =  dh/dt.  The  result  is 
referred  to  as  the  We  equation,  since  the  vertical  change  in 
h  is  due  only  to  entrainment  (no  advection) .  The  We  equation 
was  formulated  as  a  MATLAB  function,  and  solved  using  an 
intrinsic  Runge-Kutta  solution  routine,  ODE45.  The  MATLAB 
function  is  restricted  to  continuous  functions.  This 
restricts  the  model  to  deepening  scenarios. 

De  Szoeke  and  Rhines  simplified  equation  (1)  by  analyzing 
the  time  evolution  of  the  mixed  layer  deepening  process  and 
determined  which  terms  could  be  neglected  during  various  time 
regimes.    Analytic  solutions  to  the  equation  were  then 


possible  in  those  specific  time  regimes.   Figure  2.1a  is  a 
reproduction  of  de  Szoeke  and  Rhines  results,  and  2.1b  depicts 
the  results  of  a  Runge-Kutta  solution  to  equation  (1). 

De  Szoeke  and  Rhines  showed  that  initially  shear 
production  due  to  entrainment  (term  B)  can  be  ignored,  that 
turbulence  spin  up  (term  C)  dominates,  and  that  the  mixed 
layer  deepens  at  a  rate  proportional  to  u*t.  This  condition 
holds  for  the  first  100  seconds  of  deepening.  Next,  surface 
wind  stirring  (term  D)  dominates.  Deepening  continues  at  a 
rate  proportional  to  t1/3.  Term  C  begins  to  build  in 
importance  between  20  and  120  minutes,  with  the  deepening  rate 
proportional  to  t1/2.  Coriolis  (term  B)  becomes  a  factor  as  the 
mean  flow  direction  is  turned  away  from  the  wind  direction. 
After  twelve  hours,  term  B  is  no  longer  dominant.  Deepening 
returns  to  a  wind  driven  scenario  and  a  rate  proportional  to 
t1/3. 


5r- 


(  I2hrs) 


t  -sec 

Figure  2.1a   de  Szoeke  and  Rhines  solution.    Solid  curve 
ratio  D/A;  dashed  curve:  ratio  B/A;  dotted  curve:  C/A  . 
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Figure  2.1b   MATLAB  ODE45  solution,  replicating  de  Szoeke  and 
Rhines  . 


C.     IMPROVED  MODEL  PHYSICS  TO  INCLUDE  EXPLICIT  TKE 

In  the  second  model,  the  basic  We  equation  is  the  same  as 

in  the  first  model,   but  three  additional  equations  are 

included:  the  mixed  layer  temperature  equation  and  the  two 

horizontal  mixed  layer  momentum  equations.   This  allows  the 

explicit  calculation  of  the  variables  of  temperature  and 
velocity.   The  full  set  of  equations  are: 

•    °     -  e   f  3 ) 

dt      pC  h  h 

£E  =  fv+-JL- £   (4) 

dt  ph  h 


£L=-fu  +  -JL- £  (5 

dt  ph         h 


3     ccgQQh3 

2m„ut  h 

dh  = PS 

dt      (lN2h*- ut—   (1-cos  ft)   +Cnu,2h2) 
2         f2 


-W{z  =  -h) 


In  this  more  complete  model  with  explicit  momentum  and  heat 
equations,  the  ability  to  add  steady  state  upwelling  was 
included  in  the  We  equation.  This  model  is  identical  to  the 
original  de  Szoeke  and  Rhines  model  if  the  upwelling  term, 
W(z=-h),  is  set  to  zero,  and  if  the  net  surface  heat  flux(Q0) 
is  neglected.  However,  when  Q0  is  nonzero  in  the  second 
model,  Q0  directly  affects  both  the  temperature  equation  (3) 
and  the  mixed  layer  depth  equation  (6),  and  it  indirectly 


affects  momentum  because  of  the  We  terms  in  equations  (4)  and 
(5)  .  In  the  momentum  equations,  i,:  and  \,  are  components  of 
wind  stress.  AT  is  the  change  in  temperature  across  the  mixed 
layer  interface  (Figure  2.2);  similarly,  AU  and  AV  are  the 
changes  in  the  velocity  components  across  the  interface. 
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Figure  2.2   Hypothetical  temperature  profile. 

The  second  model's  solution  for  equation  (6)  was 
identical  to  de  Szoeke  and  Rhines  when  Qc  and  W(z=-h)  were  set 
to  zero.  The  three  additional  equations  allowed  the  explicit 
calculation  of  mixed  layer  temperature  and  velocity,  but  did 
not  change  the  We  solution.  The  model  was  initialized  with 
both  horizontal  velocity  components  equal  to  zero  and  an  SST 


of  17°  C  with  a  linear  temperature  gradient  with  depth,  having 
N  =  2n/600  seconds.  A  steady  wind  stress  of  0.1  N/m2  was 
applied  in  the  ix  direction.  In  30  hours,  the  SST  decreased 
to  16.3  C  (Figure  2.3)  as  the  cooler  sub-mixed  layer  water  was 
entrained.  The  horizontal  velocities,  u  and  v,  respond  to  the 
steady  state  wind  stress  in  the  form  of  forced  inertial 
motion.  The  period  is  dependent  on  latitude  (2n/f)  but  a  90 
degree  phase  shift  between  u  and  v  is  constant. 

The  third  model  is  a  MATLAB  version  of  the  NPS  mixed 
layer  model;  it  represents  a  significant  addition  of  TKE 
physics  to  the  first  two  models.  Equations  for  diurnal 
heating  and  cooling  and  three  for  energy  transformation  were 
added,  and  the  We  equation  was  reformulated.  The  resultant 
equations  are: 

Q0  =  -Q  +Q  n  (max  (0,  cos— ^— )  )  (7) 
0    a       a  86400 

-2/7?,  u?  +  U2W   +2jh,  (E-2u2)  JE-—  (  m,  E  2  +  m^fhE)         (8) 

dt  3  e      2  3     1        5 

=  V2W   +2m,  (£-3v2)  JE  -  —  (  m,  E  2  +m^fh  E)  (9) 

dt        e     2  3    1       5 

d'w2 h                              QtghQ  2     4  /in\ 

=  -cxgh&TW +2mn(E-3w2)JE--  (m,  E  2+m^fhE)  d°) 

dt      y      e    pC      2        V    3  1  5 

p 

W   =J1? (ID 

e  v    aghAT+E 
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Figure  2.3  Solution  for  temperature  and  velocity 
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$±=W   -W(z  =  -h)  (12 

dt     e 


A  Q0(t)  equation  (7)  was  designed  to  have  a  diurnal  structure 
with  a  net  zero  heat  flux.  The  energy  transformation  terms 
(m2)  redistribute  energy  among  components,  representing 
pressure-strain  interactions.  The  variable  E  is  the  summation 
of  the  squares  of  the  turbulent  velocity  fluctuations  in  the 
x,  y,  and  z  directions,  E=u2  +  v2  +  w2  ■  In  equation  (8),  m3u*3 
is  the  near  surface  wind  shear  production  term.  The  summation 
of  the  three  ir^E372  +  1% f hE  terms  found  in  equations  (8),  (9), 
and  (10)  is  the  viscous  dissipation.  The  m3u*3  term  and  the 
summed  mxE3/2  +  m5fhE  terms  are  approximated  in  the  de  Szoeke 
and  Rhines  equation  by  the  m0u,3h  term.  The  m2(E-3u  2)  E1/2  terms 
are  pressure  redistribution  of  TKE  among  the  components,  in 
which  the  turbulence  tends  toward  isotropy  by  pressure-strain 
interactions.  The  aghATWe  term  in  equation  (10)  is  the 
buoyant  damping  of  TKE  by  entrainment  and  corresponds  to  de 
Szoeke  and  Rhines  N2h4  term.  The  aghQ0/pCp  term  in  equation 
(10)  is  the  buoyant  damping  of  TKE  by  surface  heating  (Q0)  , 
and  has  no  corresponding  term  in  equation  (1),  although  it  was 
included  in  equation  (6)  of  the  second  model.  Equation  (11) 
is  the  entrainment  rate  as  a  function  of  the  ambient  TKE. 
Equation  (12)  is  the  rate  of  change  of  the  mixed  layer  due  to 
entrainment,  We,  and  vertical  mean  advection  by  W(z=-h) . 
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Initial  runs  of  this  model  were  conducted  with  a  constant 
wind  stress  but  no  heating  or  cooling.  These  solutions  were 
similiar  to  the  results  of  models  one  and  two,  but  the  mixed 
layer  depths  were  significantly  shallower  (Figure  2.4).  This 
is  as  expected  and  is  attributed  to  energy  being  expended  for 
turbulence  spin  up  and  dissipation  in  the  layer  as  well  as 
entrainment.  Two  differences  between  the  first  and  third 
models  are  significant:  (i)  calculation  of  dissipation  and 
(ii)  role  of  inertial  oscillations. 

Term  D,  in  the  de  Szoeke  and  Rhines  model,  near-surface 
wind  shear  production  of  TKE  minus  viscous  dissipation,  is  an 
over  simplification.  The  viscous  dissipation  should  not  be 
represented  as  directly  proportional  to  u.3h.  The  dissipation 
resulting  from  the  shear  production  due  to  entrainment  is  not 
represented  at  all.  Figure  2.5  is  analogous  Figure  2.1,  and 
it  shows  the  spin  up  of  model  terms.  The  ratio  of  D3/A3  is 
near-surface  wind  shear  production  of  TKE  minus  viscous 
dissipation  over  buoyant  damping  of  TKE  by  entrainment.  The 
ratio  of  B3/A3  is  shear  production  of  TKE  by  entrainment  over 
buoyant  damping  of  TKE  by  entrainment.  The  ratio  C3/A3  is  the 
turbulence  spin-up  term  over  buoyant  damping  of  TKE  by 
entrainment.  The  most  obvious  difference  is  in  the  ratio 
D3/A3.  Total  dissipation  is  much  larger  in  the  third  model 
than  with  the  approximation  used  in  the  first  model  because 
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Figure  2.4  Comparison  of  mixed  layer  depths  between  model  two 
and  model  three. 
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Figure  2.5    Model  three  solution.   The  ratios  represent  the 
same  physical  processes  as  Figures  2.1a  and  2.1b. 
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the  first  model  estimated  viscous  dissipation  as  proportional 
to  the  wind  stress  only.  The  third  model  calculates  the 
dissipation  for  all  three  TKE  components  based  on  the 
explicitly  calculated  TKE  in  the  mixed  layer.  Between  103 
seconds  and  4X104  seconds,  there  is  more  viscous  dissipation 
than  wind  shear  production,  because  of  the  addition  of  the 
dissipation  due  to  shear  production. 

The  significance  of  term  B  decreases  as  the  mixed  layer 
deepens.  In  the  first  model,  the  effect  of  inertial  forcing 
is  negligible  after  6xl04  seconds.  In  the  third  model,  term 
B3  remains  significant  four  times  longer  than  for  model  one. 
If  the  dissipation  due  to  shear  entrainment  is  neglected,  the 
amplitude  of  term  B3  decreases,  and  term  A3  and  D3  dominate 
(Figure  2.6)  after  lxlO5  seconds. 
D.    MODEL  RESOLUTION  IN  THE  TIME  DOMAIN 

Models  one  through  three  all  use  a  Runge-Kutta  solution 
with  a  variable  time  step  to  solve  the  mixed  layer  equations. 
The  numerical  scheme  maximizes  the  time  step  to  solve  the  set 
of  equations  without  exceeding  specified  error  tolerances.  In 
mixed  layer  physics,  dh/dt  is  not  a  continuous  function  when 
mixed  layer  shallowing  occurs.  For  this  reason,  the  Runge- 
Kutta  solution  can  only  be  used  during  periods  of  mixed  layer 
deepening.  To  study  realistic  cases,  a  model  that  can  shallow 
as  well  as  deepen  was  required. 
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Figure  2.6    Model  three  solution  with  viscous  dissipation 
Coriolis  parameter  m5=0  . 
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The  NPS  mixed  layer  model  was  chosen.  The  NPS  model  is 
a  FORTRAN-based  gridded  numerical  model  with  a  fixed  time  step 
of  one  hour.  Using  a  fixed  time  step  to  simultaneously  solve 
differential  equations  introduces  errors.  To  determine  the 
adequacy  of  the  one-hour  time  step,  the  NPS  model  was  compared 
to  the  third  model  using  mixed  layer  depth  (Figure  2.7)  and 
temperature  (Figure  2.8)  as  the  basis  of  comparison.  The  NPS 
model  does  not  have  a  term  for  shear  production  of  TKE  by 
entrainment  (term  B3)  but  vertical  mixing  by  dynamic 
instability  occurs  when  the  Richardson  number  (Ri)  is  less 
than  the  critical  value,  Ri  <  0.25,  following  the  procedure  of 
Adamec  et  al.  (1981) .  These  physical  processes  can  then  be 
selectively  included  or  neglected  for  comparison  purposes. 

The  NPS  model  deepens  too  quickly  during  the  first  time 
step.  This  is  because  the  NPS  model  assumes  fully  spun-up 
steady  state  turbulence.  This  is  accomplished  by  solving  the 
steady  state  equations  (13),  (14),  and  (15),  or  by  setting 
the  left  hand  side  of  equations  (8),  (9),  and  (10)  to  zero. 

7  1 

0=V2We  +  2m2(E-3v2)  JE-—  (m1E  2+m5fhE)  (13) 

aghQ  2  |  ,  ■,  * , 

0  =  -ughATW  - +2m,  (E-  3  w2)  JE-  —  {m,  E  2  +m,  fhE)  (14> 

e        pC  2  v         3       l  5 

p 

2 
0  =2m3u?  +  U2We  +  2m2  {E-3u2)  /E-  —  (m1E  2  +m5fhE)  (15) 
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Figure  2.7  Mixed  layer  depth  comparison  between  the  NPS  mixed 
layer  model  with  Richardson  Number  deepening  turned  off  and 
model  number  three  with  Term  B3  turned  off. 
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Figure  2.8  Mixed  layer  temperature  comparison  between  the  NPS 
mixed  layer  model  and  model  number  three,  at  time  =  8  hours. 
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This  assumption  allows  instantaneous  deepening  to  occur 
without  having  to  spin  up  the  turbulence.  Depending  on  the 
rate  of  fluctuation  in  the  forcing,  this  may  not  be  a  valid 
assumption.  Figure  2.9  depicts  a  portion  of  the  C3/A3  ratio 
from  Figure  2.5.  Examination  of  the  results  depicted  in 
Figure  2.9  reveals  that  the  unsteady  term  does  not  approach 
zero  at  104  seconds  as  proposed  by  de  Szoeke  and  Rhines  but  at 
105  seconds.  At  0.36xl04  seconds,  this  assumption  results  in 
a  maximum  14  percent  error  in  TKE  generation  for  any  change  in 
forcing.  If  the  forcing  is  constant  or  slowly  changing  this 
error  is  negligible.  Figure  2.9  indicates  that  the  unsteady 
term  persists  and  oscillates  in  excess  of  3  days. 

The  temperature  profiles  differ  for  two  reasons.  The 
first  is  related  to  mixed  layer  depth.  The  NPS  mixed  layer 
model  depth  is  shallower  and  therefore  warmer.  Conversely, 
model  three  has  a  deeper  mixed  layer  depth  and  is  cooler.  The 
second  reason  is  related  to  model  thermal  structure.  Model 
three  has  a  continuous  function  for  temperature  versus  depth. 
The  NPS  model  is  a  gridded  model  and  assumes  a  step  like 
function  for  temperature.  The  difference  in  thermal  structure 
between  the  models  is  evident  in  Figure  2.8  in  the  region 
below  the  mixed  layer. 

To  further  test  the  adequacy  of  the  one  hour  time  step, 
the  NPS  model  was  run  with  a  15  minute  time  step.   The  results 
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Figure  2.9  Ratio  of  Model  three  turbulence  spin-up  over 
buoyant  damping  of  TKE  by  entrainment  (C3/A3) . 
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are  depicted  in  Figure  2.10.  After  the  second  time  step  of 
the  one  hour  model,  the  mixed  layer  depths  of  the  two  NPS 
models  are  parallel  in  time  with  an  offset  of  approximately 
0.5  m.  The  one-hour  time  step  of  the  NPS  mixed  layer  model  is 
sufficient  for  the  proposed  study.  This  conclusion  was 
reached  because  the  error  induced  by  the  relatively  large  time 
step  is  small  and  explainable.  After  initial  spin  up,  the  15 
minute  time  step  does  not  result  in  significant  resolution 
improvement.  Additionally,  the  subsequent  forcing  data  to  be 
used  has  a  sampling  period  of  three  hours.  Interpolating  to 
time  steps  less  than  an  hour  was  therefore  judged  to  be 
inappropriate . 


23 


4- 


£  6 

Q. 
CD 
Q 


8- 


10 


12 


I                          I                           I                          I                          I                          I                          1 

\\ 

\\ 

\  \ 

\  \ 

\   x 

\\  \ 

\  •  \ 

-       \    ' 

— 

\    ^ 

Yw 

Model  3  -  Term  B3 

V  A 

\ 

-  -   NPS  -  Ri  deepening  (1hr) 

\x. 

NPS  -  Ri  deepening  (1/4hr) 

- 

V         N. 

N.          X. 

X              >Ss^ 

"X       ^V_ 

"S          ^>\^ 

—                        ^-^ 

1                     1                     1                     1                     1                     1                     1 

0 


8 


Time  (hrs) 


Figure  2.10  Mixed  layer  depth  comparison  using  the  one-hour 
time  step  and  the  15  minute  time  step  NPS  mixed  layer  model, 
and  model  number  three. 
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III.  REALISTIC  SOLUTIONS 

A.    IMPORTANCE  OF  VERTICAL  RESOLUTION 

In  any  geophysical  process  model,  the  resolution  required 
is  dependent  on  the  scale  of  the  physical  phenomena  to  be 
represented.  The  annual  fluctuation  in  the  mixed  layer  depth 
at  station  Papa  is  on  the  order  of  150  m.  The  daily 
fluctuations  vary  between  order  10  m  and  order  100  m, 
depending  on  time  of  year  and  the  synoptic  forcing.  Unlike 
other  mixed  layer  models,  the  NPS  mixed  layer  model  used  here 
employs  a  floating  grid  point  assigned  to  the  mixed  layer 
depth  to  increase  the  model's  resolution.  This  feature 
becomes  more  significant  as  grid  size  increases. 

The  choice  of  vertical  grid  spacing  (Az)  is  governed  by 
several  factors.  The  two  factors  considered  in  this  thesis 
are  model  computational  speed  and  resolution.  These  two 
factors  are  not  independent  of  each  other.  Smaller  grid 
spacing  requires  increased  computer  time,  but  resolves  more 
features;  larger  grid  spacing  is  faster,  but  has  reduced 
resolution.  Model  requirements  will  determine  which  one  of 
these  two  factors  is  more  important.  In  the  past,  embedding 
a  mixed  layer  within  an  ocean  general  circulation  model  or 
coupling  a  well-resolved  upper  ocean  model  to  a  global 
atmospheric  model  has  been  prohibitive  because  of  the  number 
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of  grid  levels  required  to  represent  the  ocean.   The  following 
sections  address  the  question  of  required  vertical  grid 
spacing  and  the  accuracy  of  larger  grid  sizes. 
B.    MODEL  FORCING 

The  atmospheric  forcing  utilized  was  calculated  from 
observations  by  a  Canadian  weather  ship  in  the  vicinity  of 
50°  North  latitude  145°  West  longitude,  nominally  referred  to 
as  station  Papa.  Three  hourly  weather  observations  were 
taken,  consisting  of  wind  speed  and  direction,  air  and  sea 
surface  temperature,  dew  point,  and  fractional  cloud  cover. 
The  three-hourly  observations  were  interpolated  to  hourly 
values.  Aerodynamic  bulk  formulas  derived  by  Large  and  Pond 
(1982)  and  modified  by  Oberhuber  (1993)  were  used  to  compute 
the  momentum,  sensible,  and  latent  heat  fluxes,  correcting  for 
atmospheric  stability.  Mechanical  bathythermographs (MBT)  were 
also  taken  at  irregular  intervals  at  station  Papa,  varying 
from  less  than  an  hour  to  several  days.  The  MBT  data  was 
manually  digitized  at  5  m  intervals.  The  MBT  data  are  used 
in  this  study  to  initialize  the  mixed  layer  model  and  then  for 
subsequent  verification.  Model  runs  were  conducted  for 
various  periods  between  1961  and  1969.  In  addition  to  being 
one  of  the  longest  running  continuous  oceanic  time  series,  the 
data  set  is  ideally  suited  to  testing  a  one  dimensional  model 
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because  the  horizontal  currents  are  weak  and  advection  due  to 
upwelling  is  generally  believed  insignificant  (Tabata,  1965) . 

C.  EFFECT  OF  VERTICAL  GRID  SIZE 

The  1  m  grid  NPS  mixed  layer  model  has  a  minimum  mixed 
layer  depth  of  1  m,  or  one  Az  .  Changing  the  grid  spacing 
raised  the  issue  of  what  the  minimum  mixed  layer  depth  or 
ceiling  should  be.  When  Az  was  increased  to  10  m  and  20  m, 
a  1/4  Az  minimum  mixed  layer  depth  was  found  to  be  optimal. 
The  1/4  Az  ceiling  resulted  in  a  30%  decrease  in  mean  error  in 
temperature  over  the  one-Az  ceiling.  In  the  presence  of 
strong  heating  and  no  wind  stress,  ceilings  smaller  than  1/4 
Az  resulted  in  the  mixed  layer  temperature  becoming 
unrealistically  hot.  The  same  unrealistic  temperatures  were 
also  obtained  in  the  1  m  grid  case  by  raising  the  ceiling 
above  one  Az . 

D.  VERIFICATION  WITH  OBSERVATIONS 

1 .  Model  Predicted  SST  vs  Bucket  Temperature 
Sea  surface  temperature  is  the  only  state  variable  for 
which  model  output  and  station  Papa  observations  can  be 
continuously  compared.  To  accomplish  this  task,  the  model 
mixed  layer  temperatures  were  compared  to  the  interpolated 
bucket  temperatures.  The  mean  difference  and  standard 
deviation  were  used  for  comparison. 
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E(  SST     ,    ~SSTh     .    .) 
n    v  calc  bucket  ' 

The  mean  difference  = 


n 

The  mean  difference  calculation  was  used  to  evaluate  the 
effect  of  varying  the  minimum  allowable  mixed  layer  depth  as 
described  in  section  III  C.  The  results  presented  in  Appendix 
A  provide  the  justification  for  selecting  1/4  Az  as  the 
minimum  ceiling  for  Az's  of  10  m  and  20  m.  The  eight  annual- 
period  model  runs  commenced  in  mid-March,  near  the  time  of 
maximum  mixed  layer  depth,  and  ran  for  365  days.  The  results 
of  the  mean  differences  were  unexpected  and  raised  questions 
about  optimal  spatial  resolution.  In  every  case  but  one,  the 
error  for  model  runs  with  1  m  grid  spacing  was  greater  than 
for  both  the  10  m  and  20  m  Az  cases.  This  implies  that  the 
coarser  grid  is  more  accurate.  The  standard  deviation  showed 
the  same  trend.  The  model  runs  depicted  in  Figure  3.1  are 
characteristic  of  the  other  years.  Spectra  were  calculated 
for  1  m,  10  m  and  20  m  grid  resolution  (Figures  3.2,  3.3,  & 
3.4,  respectively).  The  annual  cycle  was  not  resolved  due  to 
the  spectral  computational  length  of  the  time  series;  however, 
the  diurnal  cycle  is  evident.  The  diurnal  cycle  was  best 
resolved  by  the  10  m  grid,  then  by  the  20  m  grid,  and  lastly 
by  the  1  m  grid. 

In  an  attempt  to  understand  why  better  model  performance 
is  achieved  with  coarse  resolution,  model  tuning  sensitivity 
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Figure  3.1  Observed  bucket  temperatures  and  model  output  for 
1966. 
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Figure  3.2  Spectra  for  1  m  grid  spacing  for  1966,   using 
tunning  coefficients  ?3  =  4.6  and  AM3  =  4.5. 


30 


Cycles/Day 


Figure  3.3  Spectra  for  10  m  grid  spacing  for  1966,  using 
tunning  coefficients  P3  =  4.6  and  AM3  =  4.5. 
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Figure  3.4  Spectra  for  20  m  grid  spacing  for  1966,  using 
tunning  coefficients  P3  =  4.6  and  AM3  =  4.5. 
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tests  were  conducted.  Tuning  coefficients  P3  and  AM3  in 
equations  (13),  (14),  and  (15)  were  varied,  and  the  model 
rerun  for  each  year.  Tuning  coefficient  P3  equals  ms/n^,  and 
regulates  the  effect  of  the  Coriolis  in  the  dissipation  term; 
larger  P3  results  in  stronger  dissipation.  Tuning  coefficient 
AM3  equals  m3  and  regulates  the  efficiency  of  shear 
production.   Larger  AM3  results  in  increased  shear  production. 

Year-long  model  runs  for  eight  years  were  made  for  each 
P3  and  AM3  pair.  The  values  ranged  for  3  to  8  for  P3  and  1  to 
9  for  AM3 .  The  standard  deviation  for  each  run  was  calculated 
and  the  results  averaged.  Figure  3.5  depicts  the  average 
error  for  the  1  m  grid  spacing  cases;  Figures  3.6  and  3.7 
depict  the  average  error  for  the  10  m  and  20  m  grid  spacing 
respectively.  A  well-defined  region  of  minimum  error  is 
evident  for  the  Az  equal  to  10  m  and  20  m  cases.  This  same 
region. is  not  as  clearly  defined  in  the  Az  equal  to  1  m  case. 

Optimal  tuning  parameters  were  chosen  which  represent  the 
(P3,  AM3)  region  of  minimum  error.  The  model  was  rerun  for 
five  different  time  frames  and  for  seven  different  tuning 
coefficient  pairs .  The  first  was  a  repeat  of  the  annual 
cycle.  The  next  four  were  10  day  runs,  one  per  season.  This 
was  done  not  only  to  elucidate  differences  in  model  accuracy 
among  the  seasons,  but  also  to  determine  if  the  1  m  Az  model 
improved  over  the  larger  grid  models  for  diurnal  and  synoptic 
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Figure  3.5   Average  standard  deviation  error  (°  C)  for  1  m 
grid  spacing. 
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Figure  3.6   Average  standard  deviation  error  (:  C)  for  10  m 
grid  spacing  with  a  ceiling  of  2.5  m. 
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Figure  3.7   Average  standard  deviation  error  (°  C)  for  20  m 
grid  spacing  with  a  ceiling  of  5  m. 
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scales.  The  standard  deviations  were  calculated  and 
tabulated.  Appendix  B  is  the  complete  listing  of  results  by 
tuning  factor  pairs. 

The  results  of  the  in-depth  tuning  study  are  both 
clarifying  and  conclusive.  The  vertical  grid  spacing  accuracy 
is  dependent  on  time  of  year  and  length  of  model  run  (see 
table  1) .  During  the  spring  and  summer  when  the  forcing  is 
dominated  by  the  diurnal  cycle,  the  larger  Az's  are  more 
accurate.  During  the  synoptically  forced  deepening  period 
of  the  fall,  the  1  m  grid  spacing  more  accurately  represents 
the  bucket  temperatures.  In  the  48  cases  where  the  grid  size 
is  a  factor  when  using  10  m  vice  1  m  Az,  an  increase  in 
standard  deviation  error  of  0.06°  C  resulted,  and  an  increase 
of  0.05°  C  resulted  in  the  35  cases  when  20  m  Az  were  used. 
The  winter  accuracy  is  generally  independent  of  grid  spacing. 
2 .  Bucket  vs  Mechanical  Bathythermograph  SST 
The  NPS  mixed  layer  model  uses  BT's  to  initialize  the 
model.  After  initialization,  surface  meteorological 
observations  are  used  to  drive  the  model.  The  surface 
observations  were  taken  at  3  hourly  intervals,  but  the  MBTs 
were  deployed  at  irregular  intervals.  For  initial  conditions, 
the  model  uses  the  BT  with  the  closest  Julian  date  to  the 
start  date  of  the  run.  In  general  the  relative  vertical 
structure  of  the  surface  layer  will  be  correct,  but  a  random 
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bias  in  the  absolute  temperature  will  be  induced. 
Additionally,  a  phase  shift  is  also  likely,  as  the  model 
starts  at  00  GMT  and  the  BT  may  or  may  not  be  at  00  GMT. 
Figure  3.8  demonstrates  the  bias  between  bucket  and  BT 
temperatures  and  the  irregular  sampling  interval  of  the  BTs. 


Az 

Time 

1  m 

10  m 

20  m 

1  m 
10  m 

1  m 

& 

20  m 

10  m 
20  m 

1  m 
10  m 
20  m 

Year 

10 

33 

11 

0 

0 

2 

0 

Spring 

0 

12 

6 

0 

1 

15 

22 

Summer 

2 

7 

31 

5 

2 

7 

2 

Fall 

37 

0 

7 

1 

4 

0 

7 

Winter 

9 

1 

0 

0 

0 

2 

44 

Table  1 .  Values  represent  the  number  of  relative  minima  in 
standard  deviation  of  the  error.  The  last  four  columns 
account  for  model  runs  were  the  standard  deviations  of  the 
error  are  equal. 
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3     4     5     6 
Time  (days) 
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Figure   3.8   Example  of  bucket   temperatures  vs  mechanical 
bathythermograph  SST  at  station  Papa. 
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IV.   SUMMARY  AND  RECOMMENDATIONS 

This  study  has  evaluated  the  response  of  the  ocean  mixed 
layer  to  transient  forcing  and  the  time  and  space  resolution 
required  to  accurately  model  upper  ocean  variability.  This 
study  began  with  replication  of  de  Szoeke  and  Rhines  (1976) 
work,  using  a  Runge-Kutta  MATLAB  solution  method.  Additional 
physics  were  added  to  their  model  including:  surface  heat 
flux,  improved  viscous  dissipation,  and  anisotropy  in  the 
three-dimensional  turbulent  kinetic  energy  (TKE)  budget.  The 
surface  heat  flux  was  added  for  completeness  but  was  not 
extensively  tested  in  evaluating  the  most  rapid  transient 
responses  with  a  less  than  diurnal  period.  The  refinement  of 
viscous  dissipation  and  TKE  budget  resulted  in  modification  of 
de  Szoeke  and  Rhines'  findings.  First,  the  unsteady  TKE  term 
is  important  for  a  minimum  of  105  seconds  rather  than  only  104 
seconds.  Next,  viscous  dissipation  should  not  be  approximated 
as  being  proportional  simply  to  u*3h.  This  parameterization 
fails  to  take  into  account  buoyancy  flux  due  to  surface  heat 
flux  and  entrainment.  The  time  scale  of  the  turbulence  also 
plays  a  role  in  dissipation,  and  is  not  previously  included. 
Finally,  shear  production  due  to  entrainment  remains 
significant  for  a  minimum  of  one  pendulum-day  instead  of  the 
one-half  pendulum-day  of  the  de  Szoeke  and  Rhines'  study. 
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Expanding  the  investigation  to  model  predictions  longer 
than  a  day,  the  NPS  FORTRAN  mixed  layer  model  was  employed. 
The  required  temporal  resolution  of  mixed  layer  models  is 
dependent  on  the  phenomena  to  be  studied.  This  study  focused 
on  the  diurnal,  synoptic,  and  annual  cycles,  for  which  the 
one-hour  time  step  of  the  NPS  FORTRAN  model  is  more  than 
adequate  to  resolve.  The  one-hour  time  step  is  well  matched 
with  the  larger  vertical  grid  spacings .  Further  research  is 
required  to  determine  whether  a  one-hour  time  step  is  too 
coarse  when  used  with  1  m  grid  spacing.  If  a  time  step 
shorter  than  one  hour  is  to  be  used,  the  unsteady  term  may 
need  to  be  added  to  the  NPS  FORTRAN  model. 

Finally,  the  question  of  spatial  resolution  was 
addressed.  The  initial  hypothesis  for  this  research  was  that 
smaller  grid  spacing  would  yield  better  results  than  larger 
grid  spacing.  The  results  reported  in  Appendix  A  and  Table  1 
were  therefore  unexpected.  With  the  exception  of  the  winter 
time  frame,  the  10  m  and  20  m  grid  spacing  generally  equaled 
or  out-performed  the  1  m  grid  spacing.  Use  of  the  20  m  grid 
spacing  during  the  winter  resulted  in  an  rms  error  of  0.05°  C, 
which  is  smaller  than  the  0.1°  C  precision  of  the  bucket 
temperatures.  Overall,  the  NPS  FORTRAN  model  skill  was 
comparable  for  1  m,  10  m,  and  even  20  m  vertical  grid  spacing. 
Therefore,  the   conclusion  of  this  thesis  is  that  use  of  10  m 
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and  20  m  vertical  grid  spacing  in  global  ocean  models  and 
coupled  air-ocean  models  to  represent  the  mixed  layer  is 
accurate  and  justified. 
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APPENDIX  A.  ANNUAL  MEAN  OF  PREDICTED  TEMPERATURE  ERROR  (DEGREES  C)  FOR 
VARIOUS  MINIMUM  ALLOWABLE  MIXED  LAYER  DEPTH 


Year 

Az 

1  Az 

0.5  Az 

0.25  Az 

1.01  m 

1961 

1 

-1.5698 

-1.5537 

-1.5257 

-1.5644 

10 

-1.0555 

-1.0425 

-1.0272 

-0.9973 

20 

-1.1308 

-1.0693 

-1.0514 

-0.9948 

1962 

1 

-0.4187 

-0.3665 

-0.3338 

-0.3959 

10 

0.1393 

0.1937 

0.2296 

0.2859 

20 

-0.0428 

0.0795 

0.1181 

0.2257 

1963 

1 

-0.2282 

-0.2014 

-0.1808 

-0.2274 

10 

0.3013 

0.3154 

0.3331 

0.3608 

20 

0.1765 

0.2451 

0.2680 

0.3266 

1964 

1 

-1.2219 

-1.1993 

-1.1838 

-1.2228 

10 

-0.7613 

-0.7451 

-0.7333 

-0.7150 

20 

-0.8379 

-0.8106 

-0.7974 

-0.7618 

1965 

1 

-1.4082 

-1.3979 

-1.3831 

-1.4124 

10 

-0.9455 

-0.9213 

-0.9069 

-0.8899 

20 

-1.0784 

-1.0427 

-1.0276 

-0.9896 

1966 

1 

-0.6312 

-0.6080 

-0.5760 

-0.6246 

10 

-0.1622 

-0.1347 

-0.1143 

-0.0815 

20 

-0.3113 

-0.2362 

-0.2094 

-0.1407 

1967 

1 

-0.5077 

-0.4958 

-0.4786 

-0.5077 

10 

-0.0159 

0.0013 

0.0163 

0.0416 

20 

-0.1604 

-0.0947 

-0.0762 

-0.0289 

1968 

1 

-0.4024 

-0.3775 

-0.3423 

-0.4074 

10 

0.0859 

0.1252 

0.1555 

0.2064 

20 

-0.0529 

0.0047 

0.0373 

0.1382 
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APPENDIX  B.  RMS  TEMPERATURE  ERROR  (DEGREES  C)  VERSUS  MODEL  AZ  AND 
TUNNING  PARAMETERS  P3  AND  AM3 


P3=3.5 

AM3=2 

Winter 

1 

10 

20 

61 

0.27 

0.27 

0.27 

62 

0.17 

0.18 

0.18 

63 

0.20 

0.20 

0.20 

64 

0.14 

0.14 

0.14 

65 

0.13 

0.13 

0.13 

66 

0.14 

0.13 

0.13 

67 

0.35 

0.42 

0.46 

68 

0.19 
Summer 

0.18 

0.18 

1 

10 

20 

61 

0.63 

0.61 

0.59 

62 

0.74 

0.71 

0.69 

63 

0.65 

0.66 

0.68 

64 

0.66 

0.58 

0.61 

65 

0.72 

0.64 

0.60 

66 

0.72 

0.72 

0.71 

67 

0.72 

0.70 

0.63 

68 

0.65 
Year 

0.63 

0.63 

1 

10 

20 

61 

1.21 

0.49 

0.51 

62 

0.43 

0.95 

0.67 

63 

1.26 

0.53 

0.59 

64 

1.15 

0.48 

0.50 

65 

1.03 

0.42 

0.43 

66 

0.54 

0.58 

0.52 

67 

1.06 

0.54 

0.64 

68 

0.98 

0.45 

0.41 

Spring 

1 

10 

20 

0.28 

0.27 

0.27 

0.21 

0.20 

0.20 

0.17 

0.17 

0.17 

0.11 

0.08 

0.09 

0.22 

0.19 

0.17 

0.10 

0.09 

0.09 

0.11 

0.09 

0.09 

0.07 

0.07 

0.07 

Fall 

1 

10 

20 

0.21 

0.22 

0.24 

0.44 

0.50 

0.54 

0.24 

0.33 

0.26 

0.33 

0.42 

0.48 

0.45 

0.54 

0.44 

0.08 

0.09 

0.10 

0.16 

0.19 

0.21 

0.17 

0.17 

0.17 
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P3=4.6 

AM3=4 

5 

Winter 

Spring 

1 

10 

20 

1 

10 

20 

61 

0.27 

0.27 

0.27 

0.27 

0.27 

0.27 

62 

0.18 

0.20 

0.19 

0.21 

0.20 

0.21 

63 

0.14 

0.14 

0.14 

0.11 

0.08 

0.09 

65 

0.14 

0.14 

0.14 

0.18 

0.18 

0.17 

66 

0.14 

0.14 

0.14 

0.10 

0.09 

0.09 

67 

0.35 

0.41 

0.47 

0.10 

0.09 

0.09 

68 

0.18 
Summer 

0.18 

0.18 

0.07 
Fall 

0.07 

0.07 

1 

10 

20 

1 

10 

20 

61 

0.63 

0.62 

0.59 

0.20 

0.22 

0.24 

62 

0.74 

0.72 

0.69 

0.42 

0.47 

0.54 

63 

0.67 

0.67 

0.69 

0.22 

0.32 

0.30 

64 

0.67 

0.59 

0.59 

0.26 

0.39 

0.48 

65 

0.72 

0.65 

0.61 

0.38 

0.48 

0.38 

66 

0.73 

0.72 

0.73 

0.09 

0.09 

0.10 

67 

0.70 

0.  69 

0.65 

0.16 

0.19 

0.21 

68 

0.65 
Year 

1 

0.64 
10 

0.64 
20 

0.17 

0.17 

0.17 

61 

1.33 

0.87 

0.91 

62 

0.53 

0.43 

0.39 

63 

1.39 

0.92 

0.99 

64 

1.24 

0.82 

0.90 

65 

1.09 

0.66 

0.76 

66 

0.68 

0.35 

0.36 

67 

1.29 

0.93 

1.03 

68 

1.09 

0.66 

0.70 
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P3=7.5 

AM3=5 

Winter 

Spring 

1 

10 

20 

1 

10 

20 

61 

0.27 

0.27 

0.27 

0.27 

0.27 

0.27 

62 

0.20 

0.20 

0.20 

0.21 

0.21 

0.21 

63 

0.20 

0.20 

0.20 

0.18 

0.17 

0.18 

64 

0.14 

0.14 

0.14 

0.10 

0.08 

0.09 

65 

0.14 

0.14 

0.14 

0.18 

0.18 

0.17 

66 

0.14 

0.14 

0.14 

0.10 

0.09 

0.09 

67 

0.47 

0.50 

0.54 

0.11 

0.09 

0.09 

68 

0.18 
Summer 

0.18 

0.18 

0.07 
Fall 

0.07 

0.07 

1 

10 

20 

1 

10 

20 

61 

0.62 

0.61 

0.58 

0.22 

0.25 

0.25 

62 

0.73 

0.72 

0.68 

0.50 

0.55 

0.59 

63 

0.66 

0.66 

0.68 

0.27 

0.37 

0.27 

64 

0.65 

0.59 

0.60 

0.38 

0.50 

0.54 

65 

0.70 

0.64 

0.60 

0.46 

0.56 

0.43 

66 

0.72 

0.71 

0.71 

0.09 

0.11 

0.11 

67 

0.69 

0.69 

0.62 

0.18 

0.21 

0.22 

68 

0.64 
Year 

1 

0.63 
10 

0.63 
20 

0.17 

0.17 

0.17 

61 

0.85 

0.51 

0.50 

62 

0.54 

0.79 

0.72 

63 

0.89 

0.59 

0.66 

64 

0.77 

0.44 

0.44 

65 

0.65 

0.40 

0.42 

66 

0.49 

0.67 

0.64 

67 

0.89 

0.65 

0.71 

68 

0.75 

0.57 

0.60 
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P3=7.5 

AM3=6 

Winter 

Spring 

1 

10 

20 

1 

10 

20 

61 

0.27 

0.27 

0.27 

0.27 

0.27 

0.27 

62 

0.20 

0.20 

0.20 

0.21 

0.21 

0.21 

63 

0.20 

0.20 

0.20 

0.18 

0.18 

0.18 

64 

0.14 

0.14 

0.14 

0.10 

0.08 

0.09 

65 

0.14 

0.  14 

0.14 

0.18 

0.18 

0.16 

66 

0.14 

0.14 

0.14 

0.10 

0.09 

0.09 

67 

0.48 

0.50 

0.54 

0.10 

0.09 

0.09 

68 

0.18 
Summer 

0.18 

0.18 

0.07 
Fall 

0.07 

0.07 

1 

10 

20 

1 

10 

20 

61 

0.62 

0.61 

0.59 

0.22 

0.25 

0.25 

62 

0.73 

0.72 

0.68 

0.50 

0.55 

0.59 

63 

0.66 

0.66 

0.68 

0.28 

0.37 

0.27 

64 

0.66 

0.59 

0.60 

0.38 

0.50 

0.54 

65 

0.69 

0.64 

0.60 

0.45 

0.55 

0.41 

66 

0.72 

0.71 

0.72 

0.09 

0.10 

0.11 

67 

0.69 

0.69 

0.62 

0.18 

0.21 

0.23 

68 

0.64 
Year 
1 

0.64 
10 

0.63 
20 

0.17 

0.17 

0.17 

61 

0.89 

0.56 

0.55 

62 

0.53 

0.72 

0.68 

63 

0.92 

0.64 

0.71 

64 

0.79 

0.50 

0.54 

65 

0.68 

0.41 

0.44 

66 

0.49 

0.60 

0.59 

67 

0.96 

0.68 

0.78 

68 

0.79 

0.60 

0.62 
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P3=7 

AM3=7 

Winter 

Spring 

1 

10 

20 

1 

10 

20 

61 

0.27 

0.27 

0.27 

0.27 

0.27 

0.27 

62 

0.20 

0.21 

0.21 

0.21 

0.21 

0.21 

63 

0.20 

0.20 

0.20 

0.18 

0.17 

0.18 

64 

0.14 

0.14 

0.14 

0.10 

0.08 

0.09 

65 

0.14 

0.14 

0.14 

0.17 

0.18 

0.17 

66 

0.14 

0.14 

0.14 

0.10 

0.09 

0.09 

67 

0.44 

0.48 

0.53 

0.10 

0.09 

0.09 

68 

0.18 
Summer 

0.18 

0.18 

0.07 
Fall 

0.07 

0.07 

1 

10 

20 

1 

10 

20 

61 

0.62 

0.61 

0.59 

0.22 

0.24 

0.25 

62 

0.73 

0.72 

0.69 

0.47 

0.53 

0.58 

63 

0.66 

0.66 

0.69 

0.26 

0.37 

0.28 

64 

0.66 

0.59 

0.59 

0.35 

0.48 

0.53 

65 

0.70 

0.65 

0.61 

0.42 

0.52 

0.39 

66 

0.73 

0.72 

0.72 

0.09 

0.10 

0.11 

67 

0.69 

0.69 

0.63 

0.18 

0.21 

0.23 

68 

C.54 
Year 
1 

0.64 

10 

0.63 
20 

0.17 

C.17 

0.17 

61 

1.01 

0.66 

0.70 

62 

0.49 

0.59 

0.57 

63 

1.01 

0.81 

0.82 

64 

0.90 

0.58 

0.71 

65 

0.78 

0.46 

0.56 

66 

0.51 

0.47 

0.45 

67 

1.05 

0.84 

0.87 

68 

0.85 

0.61 

0.67 

51 


P3=8 

AM3=4 

Winter 

Spring 

1 

10 

20 

1 

10 

20 

61 

0.27 

0.27 

0.27 

0.27 

0.27 

0.27 

62 

0.20 

0.20 

0.20 

0.21 

0.21 

0.21 

63 

0.20 

0.20 

0.20 

0.18 

0.17 

0.18 

64 

0.14 

0.14 

0.14 

0.11 

0.08 

0.09 

65 

0.14 

0.14 

0.14 

0.19 

0.18 

0.17 

66 

0.14 

0.14 

0.14 

0.09 

0.09 

0.09 

67 

0.50 

0.51 

0.54 

0.11 

0.09 

0.09 

68 

0.18 
Summer 

0.18 

0.18 

0.07 
Fall 

0.07 

0.07 

1 

10 

20 

1 

10 

20 

61 

0.62 

0.61 

0.58 

0.23 

0.25 

0.25 

62 

0.73 

0.71 

0.68 

0.52 

0.56 

0.59 

63 

0.65 

0.65 

0.67 

0.28 

0.37 

0.28 

64 

0.65 

0.59 

0.61 

0.41 

0.51 

0.54 

65 

0.69 

0.64 

0.60 

0.50 

0.58 

0.45 

66 

0.71 

0.71 

0.71 

0.10 

0.11 

0.11 

67 

0.69 

0.69 

0.62 

0.18 

0.21 

0.22 

63 

0.64 
Year 
1 

0.63 
10 

0.62 
20 

0.17 

0.17 

0.17 

61 

0.73 

0.43 

0.48 

62 

0.62 

0.95 

0.82 

63 

0.80 

0.54 

0.58 

64 

0.67 

0.38 

0.37 

65 

0.56 

0.43 

0.49 

66 

0.54 

0.82 

0.78 

67 

0.85 

0.65 

0.65 

68 

0.71 

0.70 

0.64 

52 


P3=8 

AM3=5 

Winter 

Spring 

1 

10 

20 

1 

10 

20 

61 

0.27 

0.27 

0.27 

0.27 

0.27 

0.27 

62 

0.20 

0.20 

0.20 

0.21 

0.21 

0.21 

63 

0.20 

0.20 

0.20 

0.18 

0.17 

0.18 

64 

0.14 

0.14 

0.14 

0.10 

0.08 

0.09 

65 

0.14 

0.14 

0.14 

0.19 

0.18 

0.17 

66 

0.14 

0.14 

0.14 

0.10 

0.09 

0.09 

67 

0.50 

0.51 

0.54 

0.11 

0.09 

0.09 

68 

0.18 
Summer 

0.18 

0.18 

0.07 
Fall 

0.07 

0.07 

1 

10 

20 

1 

10 

20 

61 

0.62 

0.61 

0.58 

0.22 

0.25 

0.25 

62 

0.73 

0.71 

0.68 

0.52 

0.56 

0.59 

63 

0.65 

0.66 

0.68 

0.28 

0.38 

0.28 

64 

0.65 

0.59 

0.61 

0.41 

0.51 

0.54 

65 

0.69 

0.64 

0.60 

0.48 

0.57 

0.43 

66 

0.71 

0.71 

0.71 

0.10 

0.11 

0.11 

67 

0.69 

0.69 

0.62 

0.19 

0.21 

0.22 

68 

0.64 
Year 
1 

0.63 

10 

0.62 

20 

0.17 

0.17 

0.17 

61 

0.77 

0.45 

0.47 

62 

0.59 

0.83 

0.77 

63 

0.82 

0.55 

0.61 

64 

0.68 

0.41 

0.40 

65 

0.57 

0.42 

0.45 

66 

0.53 

0.73 

0.71 

67 

0.86 

0.64 

0.67 

68 

0.73 

0.67 

0.62 

53 
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